"""Present Calls

.. helpdoc::
Calculates differential expression of genes from an eSet object

"""

"""<widgetXML>
    <name>Present Calls</name>
    <icon>readcel.png</icon>
    <tags> 
        <tag>Microarray</tag> 
    </tags>
    <summary>Calculates differential expression of genes from an eSet object</summary>
    <author>
        <authorname>Kyle R Covington</authorname>
        <authorcontact>kyle@red-r.org</authorcontact>
    </author>
</widgetXML>
"""

"""
<name>Present Absent Calls</name>
<description>Calculates present and absent calls for genes from an eSet object</description>
<tags>Array Analysis</tags>
<RFunctions>panp:pa.calls</RFunctions>
<icon>readcel.png</icon>
<priority>2010</priority>
<inputWidgets></inputWidgets>
<outputWidgets>affy_diffExp</outputWidgets>
"""
from OWRpy import *
import redRGUI
import signals
class panpCalls(OWRpy):

    def __init__(self, **kwargs):
        OWRpy.__init__(self, **kwargs)
        self.require_librarys(['panp'])
        self.data = None
        self.eset = ''
        self.panpinfo = '' #used to communicate info after session reload.
        
        self.looseCut = '0.02'
        self.tightCut = '0.01'
        self.percentA = '20'
        
        self.setRvariableNames(['PA','PAcalls','PAcalls_sum','Present','peset'])
        
        """.. rrsignals::"""
        self.inputs.addInput('neset', "Normalized Eset", signals.affy.REset, self.process)
        
        """.. rrsignals::"""
        self.outputs.addOutput('pgs', "Present Gene Signal Matrix", signals.base.RMatrix)
        
        """.. rrsignals::"""
        self.outputs.addOutput('pcf', 'PA Calls Fit', signals.base.RModelFit)
        
        
        #GUI
        box = redRGUI.base.widgetBox(self.controlArea, "Options")
        
        """.. rrgui::"""
        self.looseCut = redRGUI.base.lineEdit(box, label = "Loose Cut", orientation = "horizontal")
        self.looseCut.setText('0.02')
        
        """.. rrgui::"""
        self.tightCut = redRGUI.base.lineEdit(box, label = "Tight Cut", orientation = "horizontal")
        self.tightCut.setText('0.01')
        
        """.. rrgui::"""
        self.percentA = redRGUI.base.lineEdit(box, label = "Percent Absent", orientation = "horizontal")
        self.percentA.setText('20')
        
        """.. rrgui::"""
        self.processbutton = redRGUI.base.commitButton(self.bottomAreaRight, "Process eSet", callback = self.processEset, processOnInput=True, processOnChange=True)
        
        self.R(r'''pa.calls2<-function (object = NULL, looseCutoff = 0.02, tightCutoff = 0.01, 
    verbose = FALSE) 
{
    if (is.null(object)) {
        cat("\nUSAGE:\n\tpa.calls(object, looseCutoff=0.02, tightCutoff=0.01, verbose = FALSE)\n\nINPUTS: \n\tobject - ExpressionSet, returned from expression-generating function, \n\t   such as expresso(), rma(), mas5(), gcrma(), etc. \n\tlooseCutoff - the larger P-value cutoff\n\ttightCutoff - the smaller, more strict P-value cutoff\n\tverbose - TRUE or FALSE\n\nOUTPUTS: \n\tReturns a list of two matrices, Pcalls and Pvals:\n\tPvals - a matrix of P-values of same dimensions as exprs(input object). Each\n\t   datapoint is the P-value for the probeset at the same x,y coordinates. \n\tPcalls - a matrix of Presence (P), Marginal (M), Absent (A) indicators\n\n")
        return()
    }
    if (require(affy) == FALSE) {
        stop("\npa.calls() requires BioConductor 'affy' package.\n  Currently, it is not installed. Please install and load, then\n  rerun pa.calls()\n\n")
    }
    if (verbose) {
        cat("\nInvoking function 'pa.calls'\n")
        cat("tightCutoff is ", tightCutoff, "\nlooseCutoff is ", 
            looseCutoff, "\n")
    }
    if (!is(object, "ExpressionSet")) {
        stop("\nAborting: object must be an ExpressionSet\n\n")
    }
    chip = annotation(object)
    if (chip == "hgu133b") {
        stop("\nHG-U133B is not currently supported. Supported chip types are\n\nHG-U133A and HG-U133 Plus 2.0\n\n")
    }
    if ((chip != "hgu133a") & (chip != "hgu133atag") & (chip != 
        "hgu133plus2")) {
        stop("\nAborting: chip type must be either HG-U133A or HG-U133 Plus 2.0 \n\n")
    }
    if ((tightCutoff > 1) | (tightCutoff < 0) | (looseCutoff > 
        1) | (looseCutoff < 0)) {
        stop("\nAborting: cutoffs must be between 0.0 and 1.0\n\n")
    }
    if (tightCutoff > looseCutoff) {
        stop("\nAborting: tightCutoff must be lower than looseCutoff\n\n")
    }
    if (verbose) {
        cat("Outputs will be a list containing two matrices of same dimensions as input full dataset:\n  1. Pcalls - a matrix of P/A/M indicators: \n  2. Pvals - a matrix of P-values \n  'P': P-values <= tightCutoff ", 
            tightCutoff, "\n  'A': P-values > looseCutoff ", 
            looseCutoff, "\n  'M': P-values between ", tightCutoff, 
            " and ", looseCutoff, "\n")
    }
    if ((chip == "hgu133a") | (chip == "hgu133atag")) {
        data(NSMPnames.hgu133a)
        NSMPnames <- NSMPnames.hgu133a
        rm(NSMPnames.hgu133a, envir = globalenv())
    }
    else {
        data(NSMPnames.hgu133plus2)
        NSMPnames <- NSMPnames.hgu133plus2
        rm(NSMPnames.hgu133plus2, envir = globalenv())
    }
    AllExprs <- exprs(object)
    NegExprs <- AllExprs[NSMPnames, ]
    AllExprs <- as.matrix(AllExprs)
    NegExprs <- as.matrix(NegExprs)
    cutoff_fcn <- function(x) {
        if (is.na(x))
            return("A")
        if (x <= tightCutoff) 
            return("P")
        else if (x > looseCutoff) 
            return("A")
        else return("M")
    }
    len <- length(colnames(AllExprs))
    cat("\nProcessing", len, "chips: ")
    flush.console()
    for (chip in 1:len) {
        xNeg <- NegExprs[, chip]
        xNeg <- sort(xNeg, decreasing = TRUE)
        yNeg <- seq(0, 1, 1/(length(xNeg) - 1))
        interp <- approxfun(xNeg, yNeg, yleft = 1, yright = 0)
        PV <- sapply(AllExprs[, chip], interp)
        PC <- sapply(PV, cutoff_fcn)
        if (chip == 1) {
            Pvals <- PV
            Pcalls <- PC
        }
        else {
            Pvals <- cbind(Pvals, PV)
            Pcalls <- cbind(Pcalls, PC)
        }
        cat("#")
        flush.console()
    }
    Pvals <- as.matrix(Pvals)
    Pcalls <- as.matrix(Pcalls)
    colnames(Pvals) <- colnames(AllExprs)
    colnames(Pcalls) <- colnames(AllExprs)
    outlist <- list(Pcalls = Pcalls, Pvals = Pvals)
    cat("\nProcessing complete.\n\n")
    flush.console()
    myX <- NegExprs
    maxLen <- 0
    for (i in 1:len) {
        stringLen <- nchar(colnames(AllExprs)[i])
        if (stringLen > maxLen) 
            maxLen <- stringLen
    }
    maxLen <- maxLen - 6
    if (maxLen < 2) {
        maxLen = 2
    }
    for (i in 1:len) {
        myY <- seq(0, 1, 1/(length(myX[, i]) - 1))
        myX[, i] <- sort(myX[, i], decreasing = TRUE)
        revInterp <- approxfun(myY, myX[, i], yleft = 1, yright = 0)
        revTight <- revInterp(tightCutoff)
        revLoose <- revInterp(looseCutoff)
        if (i == 1) {
            cat("\nIntensities at cutoff P-values of ", looseCutoff, 
                " and ", tightCutoff, ":\n")
            cat("Array:")
            for (j in 1:maxLen) {
                cat(" ")
            }
            cat("value at", looseCutoff, "   value at", tightCutoff, 
                "\n")
        }
        cat(colnames(AllExprs)[i], "\t", format.pval(revLoose, 
            digits = 3), "\t\t", format.pval(revTight, digits = 3), 
            "\n")
    }
    cat("\n")
    cat("[NOTE: 'Collapsing to unique x values...' warning messages are benign.]\n\n")
    return(outlist)
}''', wantType = 'NoConversion')
    def process(self, dataset):
        print 'on procress panp'
        
        if dataset == None: 
            self.status.setText("Blank data recieved")
            self.eset = ''
        else:
            self.data = dataset
            self.eset = self.data.getData()
            self.status.setText("Data Received")
            if self.processbutton.processOnInput():
                self.processEset()
            
    def processEset(self):
        if self.eset == '': return
        self.status.setText("Processing Started!!!")
        self.R(self.Rvariables['PA'] + '<-pa.calls2('+self.eset+', looseCutoff='+str(self.looseCut.text())+', tightCutoff='+str(self.tightCut.text())+')','setRData', wantType = 'NoConversion')
        self.status.setText('PA calls have been calculated')
        
        pafit = signals.base.RModelFit(self, data = self.Rvariables['PA'])
        self.rSend('pcf', pafit)
        
        
        self.R(self.Rvariables['PAcalls'] + '<-' + self.Rvariables['PA'] + '$Pcalls == "A"','setRData', wantType = 'NoConversion')
        self.R(self.Rvariables['PAcalls_sum'] + '<-apply(' + self.Rvariables['PAcalls'] + ', 1, sum)','setRData', wantType = 'NoConversion')
        self.R(self.Rvariables['Present'] + '<- ' + self.Rvariables['PAcalls_sum'] + '/length(' + self.Rvariables['PAcalls'] + '[1,]) > '+str(self.percentA.text())+'/100','setRData', wantType = 'NoConversion')
        self.R(self.Rvariables['peset']+'<-exprs('+self.eset+')[' + self.Rvariables['Present'] + ',]','setRData', wantType = 'NoConversion')
        self.R('colnames('+self.Rvariables['peset']+') <- colnames(exprs('+self.eset+'))', wantType = 'NoConversion')
        self.panpinfo = 'Processed with loose cut off = '+unicode(self.looseCut.text())+', tight cut off ='+unicode(self.tightCut.text())+', and percent absent = '+unicode(self.percentA.text())
        self.status.setText('Processed')
        
        senddata = signals.base.RMatrix(self, data = self.Rvariables['peset'])
        self.rSend('pgs', senddata)
